




-z>eA 


IMSC-HKC TR Di97S2t 


NUMERICAL ANALYSIS 
OF NATURAL CONVECTION 
IN TWO-DIMENSIONAL SQUARE 

a^:d circular containers 

IN LOW GRAVITY 


Aigvst 1980 


CoRtrpct NASW-3281 (iRterim Report) 


Prepored for 

NASA HEADQUARTERS 
WASHINGTON, DC 20546 


S. J. Robertson 

L. W. Sprodley 

M. P. Goidsteln 

Leckheed Missiles 8 Spoce CenpoRy, Irc. 
HiRtSYille Reseorch 8 ERpiRoerlRg CoRter 
4800 Bradford Drive, Hoetsville, Al 35807 



(NASA-CU-1 bJ487) NUMEuIcAi. hwALY^jIS 

NATURAL CUKVJXJTION IN TW J-UlAhNSlONAi. SijUAtiL 

ANU CxRCULAR CONTAIN EL8 iw LOw OxjAVITY 

luteriB Report (Lockheed ttio^xie:^ ana Ipaco 

Co.) J9 p tiC AoJ/MF AJ1 CaCL oj/Jw 


j8a-3068o 


(I 

1 b834 





t 


LMSC-llREC TR D697821 


■| 


FOREWORD 



This document reports tlu' results of effort by personnel 
of Lockheed Missiles Space Company, Inc., Huntsville Re- 
search liit Engineering Center, for the National Aeronautics and 
Space Administration under Contract NASW-3^81, "Manufactur- 
ing in Space: Fluid Dynamics Numerical Analysis." The NASA 
technical director for tltis contract is Dr. Robert F. Dressier, 
Manager, Advanced Technology Program, NASA Headquarters, 
Washington, D. C. 
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ABSTRACT 


A numerical study of natural convection in circular cylinder and square 
enclosures shows that the analytic low Rayleigh number theory of previous 
investigators is valid for Rayleigh numbers up to 1000. For a Rayleigh num- 
ber of 5000, steady state values of maximum fluid velocity differ by 20 per- 
cent. This deviation between analytic theory and numerical results increases 
for higher Rayleigh numbers. In addition, the low Rayleigh number theory is 
shown to be valid for higher Rayleigh numbers for a portion of the transient 
phase before significant deviation becomes apparent. It is also shown that 
square shaped experiment configurations may be analytically approxinnated 
with good accuracy by circular cylinders of equal cross sectional area for 
the prediction of convection velocities and flow patterns at low Rayleigh 
number. 
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NOMENCLATURE 


Description 

cylinder diameter = 2 R 

gravity force 

length of square side 

radial distance 

cylinder radius 

Rayleigh number = ^ , 

i/Ot 

- K3ATL^ 
po ’ 


cylinder 

square 


temperature 

initial mid-point temperature 

temperature difference across circular cylinder 
or square enclosure 

time 

dimensionless time = i/t/R^, cylinder 

2 

i^t/(L/2) , square 

velocity 

dimensionless velocity = ^ v, cylinder 

g(3 AT d^ 

— V, square 

gpAT(L) 

rectangular coordinates (Figs. 1 and 9) 
thermal diffusivity 

volumetric coefficient of thermal expansion 
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polar angle (Fig. 1) 
kinematic viscosity 
stream function 
dimensionless stream function 


— - »|j, cylinder 
gp AT d^ 


128 1/ 
g p AT 


square 
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1. INTRODUCTION 

The NASA Materials Processing in Space program includes a number 
of experiments to be performed using the Fluids Experiment System (FES) 
facility aboard the Spacelab vehicle in orbit. These experiments will investi- 
gate various fluid phenomena of interest to materials processing in space 
applications. • Fluid convection induced by residual gravity forces under 
orbital conditions is of interest to a number of possible materials process- 
ing applications. 

This report describes numerical fluid dynamics analyses of simplified 
geometries to investigate fluid convection behavior under various conditions 
of thermal gradients and gravitational loading. These combined effects are 
indicated by the magnitude of the Rayleigh number (see Nomenclature). 

The particular configurations investigated in this study were two- 
dimensional circular cylinder and square enclosures with horizontally im- 
posed initial temperature gradients. The fluid was water in both cases. 
Batchelor (Ref. 1) investigated rectangular enclosures using non-numerical 
analytical methods for the limiting case of low Rayleigh number (Ra — ► 0) 
and steady state conditions. Weinbaum (Ref. 2) made a similar study for the 
circular cylinder. Dr. Robert F. Dressier (Ref. 3) of NASA Headquarters has 
recently developed a non-numerical transient solution for the cylinder con- 
figuration. 

The purpose of the numerical study described herein is to establish 
the range of Rayleigh numbers for which the low Rayleigh number theories 
of Batchelor, Weinbaum and Dressier may be considered valid. We used 
the Lockheed developed LOCAP (LOc kheed Convection Analyzer Program) 
and GIM (General Interpolant Method) codes to perform these numerical 
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roinputalLons. Tlic LOCAP code is a two-dimensional code for analyzing 
convection in rectangular enclosures, and the GIM code is a general purpose 
fluid mechanics code for all kinds of geometries and flow conditions. The 
LOCAP computations were performed on the NASA-MSFC Univac 1108 com- 
puter, and the GIM code on the much faster and more efficient NASA -Langley 
STAR- 100 vector computer. 
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2. CIRCULAR CYLINDER ENCLOSURE 

Numerical computations were made for natural convection within the 
circular cylinder configuration shown in Fig. 1. 



Fig. 1 - Geometry for Circular Cylinder Enclosure 


The initial temperature distribution w is based on a horizontal gradient in 
the X direction, with the boundary points held constant in timet 


T(r,o) = + (AT/2) x/R 


T(R, t) = T^ + (AT/2) CO80 


( 1 ) 
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The geometry module of the GIM code divided the circular region into 
an array of 20 x 20 area elements with 21x21 nodal points. The process 
treated the boundary as a four-sided figure, each side being a circular arc. 

By use of stretching functions, each area element became a quadrilateral, 
with curvilinear sides, and with the nodal points located at its four corner 5 
From the discussion of the computer geiierated velocity vector plots, it will 
become clear where the nodal points were placed. 

We obtained transient flowfield computations over a range of Rayleigh 
numbers from 1 to 100,000. Figure 2a shows a time history of maximum 
computed flowfield velocities for various Rayleigh numbers. Both velocity 
and time are dimensionless as described in the nomenclature. 

Up through Ra = 1000, the values of v closely follow the same curve 
^ max 

through steady state. This curve is also close to the theoretical transient of 
Dressier (Ref. 3). Some deviation from the low Ra results begins to occur at 
t = 0. 1 for Ra = 5000 and 10,000. The steady state values for these cases 
were extrapolated from the derivatives of since the numerical calcula- 

tions converged slowly. 

For higher Ra the deviations begin at earlier times. Thus, for Ra = 100,000, 

the low Ra curve is followed until t = O.Oo where v peaks and then falls to a 

max ^ 

lower steady state value. Figure 2b shows a time history of deviation from 

the low Ra values. For Ra = 5000, the deviation is less than % throughout the 
transient phase, while for Ra = 10,000 this level is reached halfway titrough the 
transient. For the higher Rayleigh numbers, 20% deviation is reached about one- 
third through the transient period. 

Figure 3 shows at steady state as a function of Ra. The computed 

results closely approach the analytic value for low Ra, showing that the theory 
is valid up to Ra = 1000. The range of values 5000 < Ra < 50,000 may be viewed 
as a "transition" region. For larger Ra, the trend in v approaches the (Ra) ^ 
dependence given by Weinbaum (Ref. 2). 
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Variation of Maximum Equilibrium Velocity with Rayleigh Number for Water 
Contained in a Horizontal Cylinder 
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The low Ra theory is also valid for higher Ra during some Initial part 
of the transient. As noted earlier, for Ra = 100,000 this initial stage last 
until t = 0.06. This is further indicated by comparison of the radial distri- 
bution of velocities for Ra = 1 and 100,000 shown in Fig. 4. Not only the 
maximum velocity, but also the entire velocity field compares well up to 
t = 0.06. Note that the computed distribution approaches the low Ra theo- 
retical steady state distribution. Note also that the early velocity distribu- 
tions arc essentially linear for a considerable distance away from the center, 
thus forming a solid body rotation inner core that gradually diminishes with 
time. 


GIM code generated velocity vector plots reveal that th? flow at low 
Raleigh number remains unicellular from initial transient through steady 
slate. The velocity vector plot shown in Fig. 5 for Ra = 1 at steady state is 
seen to be nearly perfectly circular. This plot is essentially duplicated 
throughout the transient phase for all Rayleigh numbers up through Ra = 1,000. 

A slight distortion in the circular flow pattern is seen for Ra - 10,000, and a 
distinct breakup into multiple cells is noted for Ra = 100,000. This departure 
from circular flow is also indicated by the velocity contour plots in Fig. 6. 

Note the skewness that has developed in the still unicellular contours for 
Ra r. 10,000. This skewness is apparently carried over into the multiple cell 
contours for Ra = 100,000. 

Ilcat transfer at low Rayleigh numbers is dominated by conduction, and 
the effect of increasing Rayleigh number is to increase the distortion of the 
temperature contours from the initial conduction profile. This effect is shown 
in the steady state temperature contour plots in Fig. 7. Note that for Ra = 1 
the contours are nearly straight lines, corresponding to a nearly pure conduction 
temperature field. The distortion in the temperature field at higher Rayleigh 
numbers leads to an unstable condition that results in departure from circular 
flow and, for sufficiently high Rayleigh numbers, eventual breakup into multiple 
cells. This is clearly indicated in the comparison of temperature contour and 
velocity vector plots shown in Fig. 8 for various times with Ra = 100,000. At 
t rr 0.1, the temperature field is greatly distorted from the initial pure conduction 
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field, with the interior contours nearly horizontal. The flow is still uni- 
cellular, but it has become somewhat skewed. By t = 0.2, the flow has 
broken up into two cells, and this has caused the temperature contours to 
loop past the horizontal and to further distort the temperature field. By 
t = 0.3, a third cell has developed and this has tended to shift the tempera- 
ture contours back toward a more stable configuration. 
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3. SQUARE ENCLOSURE 


The square enclosure configuration is illustrated in Fig. 9. 



Hot Side 


Numerical computations were made with the LOCAP code using a 10 x 10 array 
of square elements and the GIM code using a 20 x 20 array. The LOCAP code 
was used to get initial results because it had existing capability for considering 
buoyarcy effects. The GIM code was used later after modifications were made 
to the program to include buoyancy. As in the circular cylinder, the initial 
temperature distribution was based on a horizontal gradient in the x direction. 
Calculations were made using the LOCAP code for both adiabatic (zero normal 
temperature gradient) and constant temperature boundary conditions on the top 
and bottom surfaces. The side surface temperatures were held constant in 
both cases. The GIM code calculations were for constant temperature bound- 
aries only. 
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Shown in Fig. 10 are the computed maximum flow velocities as a func- 
tion of Rayl^'igh number after the flow has reached an equilibrium condition. 

The theoretical low Rayleigh number limiting value, based on the first order 
stream function in Ref. 1, is shown for comparison. Good agreement is noted 
between the theoretical low Rayleigh number limit and the computed results. 

The LOCAP results Indicate significant deviation from the low Rayleigh num- 
ber results starting at about Ra = 1000. The GIM code results follow the same 
general trend, but the deviation occurs at higher Rayleigh number. Note that 
the adiabatic boundary produces the greatest deviation from the low Rayleigh 
number results in the LOCAP data. 

LOCAP generated plots of velocity vectors, streamlines and absolute 
velocities are shown in Figs. 11, 12 and 13. Note the skewness that develops 
in the unicellular flow patterns as Rayleigh number increases. At Ra = 100,000, 
the flow appears to be nearly broken up into two cells. LOCAP plots of tem- 
perature contours are shown in Fig. 14. Note that the adiabatic boundary 
permits a greater distortion in the temperature contours than the constant 
temperature boundary. This probably accounts for the greater deviation 
from the low Rayleigh number results shown in Fig. 10 for the adiabatic 
boundary cases. 

GIM code generated velocity vector, absolute velocity and temperature 
contour plots are shown in Figs. 15, 16 and 17 for constant temperature bound- 
aries. The GIM code results generally confirm the LOCAP data. The GIM 
code, however, predicts a definite breakup into three convection cells at 
Ra = 100,000, whereas the LOCAP code noted only an apparent incipient 
double cell pattern at the same Rayleigh number. 

The differences between the LOCAP and GIM results may be attributed 
not only to differences in the LOCAP and GIM computational techniques, but 
also to the fact that the GIM computations were based on finer grid of nodal 
points. The finer grid for the GIM computations was made possible by the 
greater speed and efficiency of the STAR- 100 vector computer. 


I6 

LOCKHEED • HUNTSVILLE RESEARCH & ENGINEERING CENTER 


i 

I , 

i . 

i •, 


' •t 



LMSC-HREC TR D697821 






f 

i 

1 

i 

( 

j 

i 


LOCKHEED • HUNTSVILLE RESEARCH & ENGINEERING CENTER 


Fig. 10 - Variation of Maximum Equilibrium Velocity with Rayleigh Number 
for Water Contained in Square Enclosure 
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The variation in maximum velocity as a function of time is shown in 
Fig. 18 for the GIM code data. The LOCAP code did not provide realistic 
results during the transient phase for the lower Rayleigh numbers. The 
GIM code results in Fig. 18 exhibit basically the same trends noted earlier 
in Fig. 2 for the circular cylinder. Again, the theoretical low Rayleigh num- 
ber limiting value, based on the first order stream function, in Ref. 1, is 
shown for comparison. 

The velocity distribution along a horizontal bisector is shown in Fig. 19 
as a function of time for Ra = 1. The low Rayleigh number theoretical equi- 
librium distribution is shown for comparison. The traversing of the theo- 
retical equilibrium curve over the GIM code results near equilibrium for 
x/(L/2) > 0.7 is obviously not realistic. The solution for the first order 
stream function given by Batchelor in Ref. 1 is actually an approximate solu- 
tion intended to provide a simple algebraic equation form for the stream 
function. Hence, exact agreement is not expected. 

A comparison of temperature contour and velocity vector plots at 
various times is shown in Fig. 20 for Ra = 100,000. A strong correlation is 
indicated between distortion of the temperature field and departure from 
symmetrical unicellular flow, as was shown earlier in Fig. 8 for the 
circular cylinder. 

A strong similarity has been noted throughout this section between the 
square enclosure and the previously discussed circular cylinder results. 

This suggests the possibility of using the simpler, more analytically defined 
circular cylinder theory to predict convection velocities and flow patterns 
in actual experiment configurations, which are more likely to be square than 
circular. In order to match velocities in a model circular cylinder configura- 
tion to an experimental square enclosure, we set the theoretical equilibrium 
velocities for the two configurations equal and find the ratio of the circle 
diameter to square side (page 29); 
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KBATd^ gQATL^ 

.385 - = .Sn 


d/L = 1.1543 

Compare this ratio to the ratio 1.1284 that yields equal areas for the circle 
and square. This somewhat fortuitous result clearly implies that an experi- 
mental square enclosure may be closely modeled by a circular cylinder of 
equal cross sectional area. The validity of this modeling is additionally 
supported by the rather remarkable degree of radial symmetry that exists 
in the low Rayleigh number velocity field in the square enclosure. This is 
displayed in the angular variation of velocity shown in Fig. 21 for various 
radial distances. These velocities were determined from the first order 
stream function in Ref. 1. 
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4. CONCLUSIONS 

The primary conclusion to be drawn from this study is that the low 
Rayleigh number theory for convection in circular cylinder and square box 
enclosures is valid throughout the transient phase to equilibrium for Rayleigh 
numbers up to at least 1000. The GIM code results, which we assume to be 
most valid, show a deviation of 20% from the low Rayleigh number results at 
Ra = 5,000, Additionally, for higher Rayleigh numbers, the low Rayleigh 
number theory is valid for some portion of the transient, depending on the 
value of the Rayleigh number, before significant deviation becomes apparent. 

In the case of both circular cylinder and square box enclosures at Ra = 100,000 
the valid portion of the transient phase extends up to"t S' 0.06, which is approxi- 
mately 20% of the time required to reach equilibrium. 

Another conclusion drawn from the comparison of the circular cylinder 
and square box enclosure results is that experimental square shaped con- 
tainers can be analytically modeled rather closely by circular cylin.’ers of 
equal cross sectional area for the prediction of convection velocities and flow 
patterns at low Rayleigh number. 
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